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A  Theory  for  the  Scalar  Roughness  and  the 
Scalar  Transfer  Coefficients  Over  Snow  and  Sea  Ice 


EDGAR  L  ANDREAS 


INTRODUCTION 

In  the  atmospheric  surface  layer  at  neutral  stability,  the  velocity  (U),  potential  tempera¬ 
ture  (T),  and  water  vapor  density  (Q)  profiles  have  the  familiar,  semi-logarithmic  form, 


^  -  k  '  ln(z/z0). 

(1) 

u. 

T(z)-T0 

-  (aH  k)-'  ln(z/zT), 

*  * 

(2) 

QUhQo 

-  (aE  k)  '  ln(z/zQ). 

(3) 

Here  z  is  the  height  above  the  surface;  k  is  von  Karman’s  constant  (0.4);  T0  is  the  surface 
temperature;  Qo  is  the  water  vapor  density  of  air  in  saturation  with  a  snow  or  sea  ice  surface 
at  T0;  and  aH  ( =  KH/KM)  and  aE  ( =  KE/KM)  are  the  ratios  of  the  scalar  turbulent  diffusivi- 
ties,  Ah  and  KE,  to  the  turbulent  diffusivity  for  momentum  (e.g.,  Dyer  1974).  The  w„  tt 
and  q ,  relate  the  profiles  to  the  turbulent  surface  fluxes  of  momentum  (t)  and  sensible  (Hs) 
and  latent  (H\)  heat; 


T  =  QU!,, 

(4) 

=  -QC  pH.'., 

(5) 

(6) 

where  e  =  air  density 

cp  =  specific  heat  of  air  at  constant  pressure 
Ls  =  latent  heat  of  sublimation  of  ice. 

Equations  1-3  define  the  roughness  lengths,  zq  is  the  familiar  roughness  length  for  wind 
speed;  zT  and  Zq  are  the  roughness  lengths  for  temperature  and  water  vapor — the  so-called 
scalar  roughness  lengths,  z0  is  the  height  at  which  the  semi-logarithmic  velocity  profile  ex¬ 
trapolates  to  U  =  0.  Similarly,  zT  and  zQ  are  the  heights  at  which  the  semi-logarithmic  tem¬ 
perature  and  water  vapor  profiles  extrapolate  to  the  surface  values,  Tq  and  Qo,  respectively. 
All  are  fictitious  levels  since  the  semi-logarithmic  profiles  are  not  valid  clear  down  to  the 
roughness  lengths. 


Knowing  the  roughness  lengths  is  equivalent  to  knowing  the  bulk-aerodynamic  transfer 
coefficients  for  momentum  (CD,  the  drag  coefficient)  and  for  the  scalars,  sensible  (CH)  and 
latent  (CE)  heat.  After  specifying  a  reference  height  r  (henceforth  taken  as  10  m)  to  be  the 
level  where  average  values  of  wind  speed  (t/r),  temperature  (7"r),  and  humidity  (Qr)  are  meas¬ 
ured,  we  define  these  transfer  coefficients  as 


r  =  qCdUt\ 

H%  =  Q ^pC^UT  (Tq-Tt), 

Hl  =  LsCEUT(Qo~Qr). 

For  neutral  stability  eq  1-6,  in  turn,  relate  these  coefficients  to  the  roughness  lengths: 
k ! 


Cd  - 

C\\  — 

Ci  = 


[\n(r/zo)Y  ' 
«h*Cd 

£Q^-ln(zT/z0) 

«E*Cg 

k  Cqa  -  ln(zQ/z0) 


(7) 

(8) 
(9) 

(10) 

(11) 

(12) 


Since  correcting  the  transfer  coefficients  for  stability  effects  is  straightforward  (e.g.,  Dear- 
dorff  1968,  Large  and  Pond  1982),  from  here  on  all  of  my  references  to  transfer  coefficients 
will  be  to  these  neutral-stability  ones.  Clearly,  from  eq  10-12,  CH  =  CD  only  when  Zj  =  z<> 
and  aH  =  1;  and  CE  =  CD  only  when  zQ  =  z0  and  aE  =  1.1  will  show  shortly  that,  contrary 
to  the  common  assumption  (e.g.,  Paulson  1970,  Businger  et  al.  1971 ,  Lettau  1979),  zT  and  zq 
rarely  equal  z0- 

A  typical  goal  of  micrometerology  is  to  find  CD,  CH  and  CE  or,  equivalently,  z0,  zT  and  zQ. 
These  are  fairly  well  known  over  the  ocean  but  are  still  only  poorly  known  over  most  other 
horizontally  homogeneous  surfaces.  In  particular,  only  CD  is  well  known  over  sea  ice.  The 
ex»-  sive  set  of  CD  values  that  Banke  et  al,  (1980)  reported  show  that  over  sea  ice  CDis  an  in¬ 
creasing  function  of  the  surface  roughness.  The  roughness  parameter  here  is  the  root-mean- 
square  (rms)  surface  elevation  along  a  line  parallel  to  the  wind  direction  and  should  not  be 
confused  with  the  roughness  length  zo-  Leavitt  et  al.  (1977)  also  found  that  CD  increased  as 
the  sea  ice  surface  got  rougher.  Arya  (1973,  1975)  had  theoretically  predicted  this  increase  in 
CD  with  roughness,  showing  it  to  be  a  consequence  of  the  form  drag. 

I  know,  however,  of  only  two  published  attempts  to  measure  CH  and  CE  over  snow  or  sea 
ice,  those  by  Hicks  and  Martin  (1972)  over  snow-covered  Lake  Mendota  and  by  Thorpe  et  al. 
(1973)  in  the  Beaufort  Sea  and  in  Robesi  n  Channel.  And  the  results  are  inconclusive.  Hicks 
and  Martin  (1972)  found  CH  =  0.9  x  10"3  and  CE  =  2.5  x  10  ',  while  Thorpe  et  al.  (1973)  re¬ 
ported  CH  =  1.2x10"’  and  CE  =  0.55x10  ’.  That  is,  in  one  case  CH/Ct  =  0.4,  and  in  the 
other  CH/CE  =  2.2.  With  only  these  few,  contradictory  values  and  without  an  adequate 
theory  from  which  to  estimate  CH  and  CE  in  the  absence  of  experimental  values,  sea  ice  mod¬ 
elers  have  had  to  rely  on  intuition  or  convention.  Most  (e.g.,  Parkinson  and  Washington 
1979,  Hibler  1980)  have  followed  Maykut  (1978)  and  assumed  that  CM  and  C,  are  constant  — 
both  equal  to  1.75  xlO'3.  But  eq  11  and  12  imply  that  CH  and  Cj  are  not  constant;  they  de¬ 
pend  on  the  characteristics  of  the  surface  and  on  the  wind  speed,  since  C|}  depends  on  the 
surface  characteristics.  Andreas  and  Ackley  (1982)  evidently  were  the  first  to  point  this  out. 


In  this  paper  I  will  present  a  theoretical  model  for  predicting  CH  and  Ct  over  snow  and  sea 
ice  that  relies  on  the  empirical  dependence  of  CD  on  surface  roughness  reported  by  Banke  et 
al.  (1980).  From  eq  1 1  and  12  it  is  clear  that  to  predict  CH  and  Ct  I  must  also  find  Zj/zq  and 
Zq/zo ■  To  do  this  I  derive  the  scalar  profiles  in  the  interfacial  sublayer  over  both  aerodynam- 
ically  rough  and  aerodynamically  smooth  surfaces.  Matching  these  to  the  semi-logarithmic 
or  inertial  sublayer  (Tennekes  and  Lumley  1972,  p.  147)  profiles  (eq  2  and  3),  I  then  solve  for 
the  scalar  roughness.  I  treat  snow  and  sea  ice  with  the  same  model  because  sea  ice  is  generally 
snow  covered;  the  two  surfaces  are  therefore  aerodynamically  similar. 


AERODYNAMICALLY  ROUGH  SURFACE 

With  the  ratio  of  kinematic  viscosity  to  molecular  diffusivity  (o),  the  roughness  Reynolds 
number  7?.  =  u,zq/c  customarily  parameterizes  wall-bounded  shear  flows.  Although  two 
flows  may  have  different  velocity  or  length  scales  or  different  kinematic  viscosities  (r),  they 
are  dynamically  similar  if  their  roughness  Reynolds  numbers  and  their  a  values  are  the  same. 
Three  dynamic  regimes  are  possible,  each  characterized  by  a  different  /?,  range  (e.g.,  Bus- 
inger  1973).  If  /?,  <  e2  =  0. 135,  the  surface  roughness  elements  are  imbedded  in  the  viscous 
sublayer  and  the  surface  is  aerodynamically  smooth.  If  /?,  >  2.5,  the  roughness  elements 
poke  through  the  viscous  sublayer  and  the  surface  is  aerodynamically  rough.  For  0.1 35  <  /?, 
<  2.5,  the  surface  is  in  transition. 

Brutsaert  (1975a)  and  Liu  et  al.  (1979)  based  models  of  the  turbulent  transfer  over  aerody¬ 
namically  rough  surfaces  on  a  surface-renewal  model  (Danckwerts  1970,  p.  100).  Small  ed¬ 
dies  continually  sweep  into  the  interfacial  sublayer,  remain  in  contact  with  the  surface  for  a 
short  time,  transferring  heat  and  moisture  by  molecular  diffusion,  and  then  finally  burst  up¬ 
ward  ahead  of  inrushing  new  eddies.  Grass  (1971)  suggested  that  while  the  eddies  are  in  con¬ 
tact  with  the  surface,  they  may  be  stagnant — trapped  by  the  roughness  elements.  Brutsaert 
(1975a)  thus  assumed  that  over  a  rough  surface  the  interfacial  transfer  of  scalar  properties  is 
strictly  a  diffusion  process.  That  is,  using  temperature  as  an  example, 


•II 

(it 


=  D 


d!T 

dz: 


(13) 


where  t  is  time.  In  eq  13  and  in  all  that  follow,  we  could  use  water  vapor  or  any  other  scalar 
that  obeys  the  same  conservation  equation  as  temperature  (Hill  1978);  the  only  changes 
would  be  in  the  molecular  diffusivity  D  and  in  the  other  thermodynamic  constants,  such  as 
the  or ’s  or  cp  and  Fv.  The  boundary  conditions  on  eq  13  are 

F  =  7“b  for  z  >  0,  t  =  0, 

T  =  7b  for  large  z,  t  >  0,  (14) 

T  =  7"0  for  z  =  0,  t  >  0, 

where  Fb  is  the  “bulk”  temperature  above  the  interfacial  sublayer.  Many  standard  texts 
show  how  to  solve  eq  13  with  the  boundary  conditions  of  eq  14  (e.g..  Duff  and  Naylor  1966, 
p.  1 18);  the  solution  is 


T(z,t)  =  (F0-Fb)  erfc 


2  (Dry 


+  Fb, 


(15) 


where  erfc  is  one  minus  the  error  function,  erf  (Abramowitz  and  Stegun  1965,  p.  297). 
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Equation  15  models  the  diffusion  into  a  single  eddy.  Because  eddies  continually  sweep 
over  the  surface,  we  must  integrate  in  time  to  find  the  average  interfacial  sublayer  tempera¬ 
ture  profile.  Brutsaert  (1975a)  and  Liu  and  Businger  (1975)  used  Danckwerts’  (1951)  distri¬ 
bution  function, 

MO  =  t}'  exp (~t/tT),  t  >  0,  (16) 

to  model  the  fraction  of  surface  area  that  has  had  eddies  in  contact  for  a  time  t.  Here  tr  is  a 
time  scale  yet  to  be  specified.  The  time-averaged  temperature  profile  is  thus 


T(z)  =  j  T(z,t)  MO  dt.  (17) 

o 

Abramowitz  and  Stegun  (1965,  p.  303)  show  how  to  integrate  the  error  function  in  eq  15;  the 
solution  of  eq  17  is  thus 

T(Z)  =  7b-(7b-rb)(l-exp(-z/6T)],  (18) 

where 

6t  =  (Dtc)'A.  (19) 


Khundzhua  and  Andreyev  (1974)  verified  eq  18  experimentally  in  the  aqueous  sublayer  in 
the  Black  Sea  and  related  the  length  scale  6T  to  the  sensible  heat  flux.  For  air  this  relation  is 


6t  =  ecpD(T0-Tb)/Hs.  (20) 

Actually,  eq  20  is  a  necessary  consequence  of  eq  18.  Hs  is  related  to  the  temperature  gradient 
evaluated  at  the  surface  by 


-QCpD 


dT 

dz 


z  = 0  ' 


(21) 


But  from  eq  18,  dT/dz\z=o 's  simply  (7b-7o)/6T;  eq  20  thus  follows  from  eq  21. 
Substituting  eq  5  for  Hs  in  eq  20,  we  can  write 


T<y-Tb  =  -u,ttdr/D. 


(22) 


Substituting  this  into  eq  18  yields  a  form  for  the  temperature  profile  in  the  interfacial  sub¬ 
layer  that  is  compatible  with  the  inertial  sublayer  profile. 


T(z)-r0 

K 


(K,  bj/D)  [l-exp(-f)]. 


(23) 


Here  f  =  z/6T. 

Brutsaert  (1975a)  and  Liu  et  al.  (1979)  assumed  that  the  time  scale  tr  in  6T,  eq  19,  is  propor¬ 
tional  to  the  Kolmogorov  time  scale  where  t  is  the  dissipation  rate  of  turbulent  kin¬ 

etic  energy.  Formalizing  this  assumption  by  defining  a  proportionality  constant  G  and  set¬ 
ting  t  =  u\/zo  (Liu  et  al.  1979),  we  have 


& 


I 
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K 


tT  =  G'{zo  v/ui)y\ 


(24) 


or 


&T/Z0  =  G  Pr  !  /? 


-l/l  O- J/i 


(25) 


where  Pr  =  v/D  is  the  Prandtl  number.  Liu  et  al.  (1979)  estimated  C  =  9.3  on  the  basis  of 
data  reported  by  Mangareila  et  al.  (1973)  for  flow  over  wind  waves.  Since  snow  or  sea  ice 
surfaces  are  less  compliant  than  water,  this  value  of  G  may  not  be  appropriate.  In  fact,  in 
evaluating  G,  Liu  et  al.  (1979)  also  considered  data  reported  by  Chamberlain  (1968)  that  im¬ 
ply  G  =  5.6  for  flows  over  various  solid  surfaces.  My  model,  with  this  G  value,  fits  data  col¬ 
lected  over  various  solid  surfaces  much  better  than  with  G  =  9.3.  Evidently,  the  proportion¬ 
ality  constant  in  eq  24  and  25  depends  on  whether  the  surface  is  firm  or  compliant. 

With  eq  25  we  can  now  simultaneously  solve  eq  2  and  23  to  find  zT.  This  is  where  1  diverge 
from  Brutsaert  (1975a).  I  will  simply  match  both  the  temperature  profiles  and  their  first  de¬ 
rivatives  at  z  =  h.  The  temperature  and  the  heat  flux  will  therefore  both  be  continuous  from 
the  interfacial  to  the  inertial  sublayer.  With  these  profile  and  first-derivative  equations,  1  can 
find  the  two  unknowns  h  and  Zj.  Brutsaert  (1975a)  never  computed  the  interfacial  sublayer 
profile;  instead  he  solved  for  zT  by  matching  the  interfacial  and  inertial  sublayer  fluxes  at  the 
arbitrary  level  h  =  7.39zo.  Although  they  do  not  say  explicitly,  Liu  et  al.  (1979)  used  the 
method  of  solution  that  I  am  proposing.  They,  however,  set  G  =  9.3  and  aH  =  aE  =  1.14, 
while  I  use  G  =  5.6  and  aH  =  aE  =  1.0. 

Matching  the  profiles  eq  2  and  23  at  z  =  h  gives 


In  Mn(zo/5T)-ln(zT/Zo)  =  A[l-exp(-f)J, 
where 


(26) 


r  =  h/h  T, 


(27) 


and 


A  =  aHkGPtVlRY\ 

A 

Matching  the  first  derivatives  at  f  yields 
fexp(-f)  =  A-\ 


(28) 


(29) 


I  do  not  need  to  worry  here  or  in  eq  26  about  the  stability  of  the  atmospheric  surface  layer 
and  its  effects  on  the  inertial  sublayer  profiles.  We  will  see  that  the  matching  level  is  well  be¬ 
low  the  region  where  atmospheric  stability  affects  the  semi-logarithmic  profiles  (Bradley 
1972).  Notice,  eq  29  has  a  solution  only  for  A  >  2.72.  Substituting  it  into  eq  26  gives  a  for¬ 
mal  expression  for  z^/zq. 


z  t/2o  =  f(6T/zo)exp(r'-A). 


(30) 


Here  A  and  Ar  are  functions  only  of  Pr  and  R ..  We  solve  for  f  by  using  Newton’s  method 
to  find  the  zeroes  in  eq  29  as  functions  of  Pr  and  R,.  Figure  1  shows  h/zo  =  (s  <5T/io)  for 
both  temperature  and  water  vapor  when  R,  =  10.  The  values  are  much  smaller  than  the  con¬ 
stant  that  Brutsaert  (1975a)  used.  He  chose  h  =  7.39z0  because  this  is  approximately  the 
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Figure  1.  Matching  height 
over  an  aerodynamically 
rough  surface  as  a  function 
of  roughness  Reynolds  num¬ 
ber.  For  temperature ,  a  = 
0. 71;  for  water  vapor,  a  =  0.63. 


Figure  2.  Matching  of  interfacial  and  inertial 
sublayer  profiles  of  temperature  (o  =  0. 71) 
and  water  vapor  (o  =  0.63)  over  an  aerody¬ 
namically  rough  surface  for  R .  =  10. 


height  of  the  roughness  elements.  By  implying  that  the  roughness  elements  protrude  above 
the  interfacial  sublayer.  Figure  1  is,  thus,  consistent  with  our  conceptual  model  of  molecular 
diffusion  into  stagnant  eddies. 

Figure  2  shows  the  transition  of  temperature  and  water  vapor  profiles  from  the  interfacial 
to  the  inertial  sublayer.  Both  the  profiles  and  their  first  derivatives  are  continuous.  The  verti¬ 
cal  fluxes  of  sensible  and  latent  heat  are,  therefore,  also  continuous. 

Figure  3  compares  predictions  of  my  model  for  the  scalar  roughness,  zs,  with  experimental 
data  from  Owen  and  Thomson  (1963)  and  Chamberlain  (1966,  1968).  All  the  data  sets  are 
from  wind  tunnel  measurements  and  represent  three  different  a  values.  Chamberlain  (1966) 
studied  water  vapor  transfer  over  toweling  and  artificial  grass  ( o  =  0.62  at  -  20°C).  Because 
the  grass,  however,  had  a  roughness  length  of  1 .0  cm— much  larger  than  that  typical  of  snow 
or  sea  ice  (Untersteiner  and  Badgley  1965,  Banke  et  al.  1980,  Schmidt  1982)— and  roughness 
elements  unlike  those  of  snow  or  sea  ice,  I  have  not  included  those  data.  Chamberlain  (1968) 
collected  his  thorium-B  ( o  =  2.78)  and  water  vapor  data  over  a  host  of  two-  and  three- 
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Figure  3.  Model  predictions  for  an  aerodynamically  rough 
surface  compared  with  measured  scalar  roughness  lengths  for 
water  vapor,  thorium-B  and  camphor. 


dimensional  roughness  elements,  some  with  roughness  lengths  as  large  as  0.6  cm.  In  the  fig¬ 
ure  I  have  indicated  the  data  for  surfaces  with  zo  >  0.2  cm — roughly  the  maximum  sea  ice 
roughness  that  Banke  et  al.  (1980)  reported.  No  systematic  difference  between  large  and 
small  roughness  lengths  is  evident,  however.  Owen  and  Thomson  (1963)  looked  at  the  trans¬ 
fer  of  camphor  (o  =  3.2)  over  glass  surfaces  with  two-  and  three-dimensional  roughness. 

Measuring  zs/zo  is  difficult  under  any  circumstances — even  in  a  wind  tunnel — because,  as 
eq  1  and  2  show  with  temperature  for  example,  we  have  to  know  u,,  t,  and  zo-  Chamberlain 
(1968)  explained  that  his  zo  values  alone  may  have  been  in  error  by  50%.  The  scatter  of  the 
data  in  Figure  3  is,  thus,  not  surprising.  In  view  of  this  uncertainty,  the  model  predictions  are 
quite  good.  The  model  reproduces  the  R, -dependence  at  constant  a  very  well  and  has  zs/zo 
decreasing  with  increasing  a,  as  the  data  do.  Only  the  zs/zo  data  from  Owen  and  Thomson 
(1963)  deviate  significantly  from  model  predictions.  1  can  only  speculate  that  significant 
measurement  errors  crept  into  their  results.  For  example,  Owen  and  Thomson  (1963)  deter¬ 
mined  the  camphor  flux  by  weighing  a  test  section  of  the  floor  of  their  wind  tunnel  before 
and  after  each  experimental  run.  The  7-cm-square  glass  test  section  must  have  been  at  least 
100  g,  but  the  typical  difference  in  camphor  coating  before  and  after  each  of  their  runs  was 
only  0.1  g.  Their  flux  measurement,  therefore,  likely  had  low  precision.  Secondly,  as  cam¬ 
phor  was  subliming  from  the  floor  of  the  wind  tunnel  during  a  run,  the  zo  value  may  have 
been  changing. 
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Figure  4.  Model  predictions  for  an  aerodynamical- 
ly  rough  surface  compared  with  the  experimental 
data  of  Dipprey  and  Sabersky  (1963). 


Figure  4  compares  model  predictions  of  z%/zq  with  the  data  of  Dipprey  and  Sabersky 
(1963),  who  investigated  heat  transfer  in  water-filled  pipes.  Pipe  flow  may  not  seem  to  be  a 
good  experimental  model  for  flow  in  the  atmospheric  surface  layer,  but  the  two  are,  in  fact, 
mathematically  equivalent.  Flow  in  pipes  is  characterized  by  a  viscous  sublayer  near  the  wall 
and  a  semi-logarithmic  inertial  sublayer  further  from  the  wall  (Schlichting  1968,  p.  578;  Ten- 
nekes  and  Lumley  1972,  p.  149),  just  like  my  model  for  the  atmospheric  surface  layer.  Figure 
4  confirms  the  validity  of  the  comparison.  Dipprey  and  Sabersky  (1963)  varied  a  by  changing 
the  water  temperature;  the  model  fits  their  data  extremely  well  for  a  =  1.20  and  reasonably 
well  at  the  other  a  values.  For  all  of  the  a  values  in  the  figure,  the  difference  between  the  data 
and  the  model  predictions  tends  to  decrease  as  R ,  increases.  This  suggests  some  experimental 
imprecision  at  low  flow  rates.  Nevertheless,  as  in  Figure  3,  the  model  does  reproduce  the  gen¬ 
eral  trends  in  the  data.  In  summary,  my  model  seems  to  be  an  adequate  fit  to  the  available 
data  that  are  most  representative  of  an  aerodynamically  rough  snow  or  sea  ice  surface. 


AERODYNAMICALLY  SMOOTH  SURFACE 

Brutsaert  (1975a)  also  modeled  scalar  transfer  over  an  aerodynamically  smooth  surface  by 
again  postulating  a  surface-renewal  mechanism.  Over  a  smooth  surface,  however,  an  im¬ 
pinging  eddy  remains  in  motion;  the  transfer  is  thus  governed  by  an  advective  diffusion 
model, 


«•  rl 


(31) 
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where  i/(z)  is  the  velocity  profile  in  the  viscous  sublayer.  The  boundary  conditions  on  eq  31 
are 


Kf 

II 

for  z  >  0, 

x  =  0, 

“i 

11 

for  large  z, 

x  >  0, 

(32) 

11 

t- 

for  z  =  0, 

x  >  0. 

In  the  viscous  sublayer,  z  u,/v 

<  5, 

V(,z)  =  «i  z/v  (33) 

(e.g..  Monin  and  Yaglom  1971,  p.  270;  Brutsaert  1975a).  On  substituting  eq  33  into  eq  31 
and  making  the  change  of  variables 


uj  zl 
V  ~  9vD  x 

suggested  by  Kestin  and  Persen  (1962),  we  get  the  equation 


,  ...  dT  d‘T 

(v+  /3)  at  +  v  -W  =°- 

The  solution  that  satisfies  the  boundary  conditions  is  (Kestin  and  Persen  1962) 

T(v)  =  7b-(7b-7b)  , 


(34) 


(35) 


(36) 


where  T  is  the  gamma  function  and  7  is  the  incomplete  gamma  function  (Abramowitz  and 
Stegun  1965,  p.  255  and  260,  respectively).  Although  Brutsaert  (1975a)  posed  eq  31  with  the 
boundary  conditions  of  eq  32,  he  solved  only  for  dT/dz  at  z  =  0. 

We  next  impose  the  assumptions  of  the  surface-renewal  model— that  the  eddy  is  in  contact 
with  the  surface  only  for  0  <  x  <  xq.  Averaging  T(rj)  over  this  distance  yields 


T(z,x 0)  =  7o- 


(7b- 7b) 

rc/3) 


[7(  Vi ,  7o)  +  70  r<-  2/j  ,  i?o)J . 


(37) 


where 


=  r(-2/j)-7(-2/3,j?0) 

is  another  form  of  the  incomplete  gamma  function,  and 


(38) 


Vo  = 


u\  zl 


9 1 D  x0 


(39) 


To  remove  the  explicit  dependence  on  x0  >n  ecl  37,  we  have  to  average  over  all  possible 
values  of  xq.  Brutsaert  (1975a)  suggested  setting 


*0  =  a  u,  /, 


£v 


(40) 


where  a  is  a  constant,  and  then  using  eq  16  for  the  distribution  function  of  /.  In  other  words, 
we  would  again  have  the  integral  eq  17,  with  eq  37  substituted  for  eq  15.  If  Brutsaert  had  at¬ 
tempted  this  integration,  he  would  have  found  the  integral  infinite.  The  reason  is  that  the  dis¬ 
tribution  function  (eq  16)  is  not  the  appropriate  one  over  a  smooth  surface.  The  work  of  Kim 
et  al.  (1971)  suggested  that  for  smooth  surfaces  the  eddy-contact  time  has  the  distribution 

<M0  =  (f/fj)  exp(-//ts),  t  >  0,  (41) 

where  rs  is  a  new  time  scale.  The  average  contact  time  is  thus  7  =  2/s.  In  addition,  from  Fig¬ 
ures  23  and  24  in  Kim  et  al.  (1971)  I  derived 

7  =  170  v/u2„  (42) 


or 


/s  =  85  v/u\.  (43) 

Notice,  since  u/u.  is  the  appropriate  scaling  length  in  the  viscous  sublayer  over  a  smooth  sur¬ 
face  (Tennekes  and  Lumley  1972,  p.  152),  v/u\  is  the  only  reasonable  time  scale  there. 

Substituting  eq  37  and  41  into  the  time-averaging  integral,  eq  17,  rearranging  arguments  a 
little,  and  defining 


7s  = 


9avDts 


we  derive  an  expression  for  the  average  profile  in  the  interfacial  sublayer, 


(44) 


(7b-7b)»»I  f 

T(z)  =  Tq - rTvfP  J  {7o'7(‘/3,7o)  +  r(-Vj,>?o)]7o  CXPf-ls/loMo-  (45) 

The  method  of  steepest  descent  (e.g.,  Dennery  and  Krzywicki  1967)  is  useful  for  approximat¬ 
ing  such  difficult  integrals;  it  yields 

T(z)  =  T0-  |0.9607(,/3.t)s/2) 

+  USS^C/j.^-re/j)  +  v*  exp(-7)s)]|.  (46) 

This,  to  my  knowledge,  is  the  first  derivation  of  the  interfacial  sublayer  profile  of  a  scalar 
over  an  aerodynamically  smooth  surface  that  is  based  on  surface-renewal  concepts. 

As  with  the  rough-surface  case,  we  have  to  eliminate  7b  in  eq  46  to  match  the  interfacial 
and  inertial  sublayer  profiles.  Again,  we  know  how  the  surface  flux  is  related  to  the  profile — 
simply  by  eq  21.  Using  this  and  substituting  eq  43  for  /s  in  eq  46  we  find 

Tq-  Tb  =  1 . 5 1 9(85a) 1/1  Pr  *  t, .  (47) 

Hence, 

T(z )  =  T0  +  1  ^^Sa)"  Pr*  r,  K(Vs),  (48) 


where 


*(t>s)  =  re/,)-'  |7('/j,v2)  +  i.44he/,,r,s)-re/,)  +  ^  e\P(-v,m. 


(49) 


Notice,  with  /s,  eq  43,  substituted  into  vs.  eq  44,  and  recognizing  that  over  a  smooth  surface 


zo  =  e v/u, 

(Tennekes  and  Lumley  1972,  p.  157),  becomes 

Matching  the  profiles  at  z  =  h  or  at 

*■  -  Pr  (|)’  • 

we  get 

'/,  In  +  ‘/s  -  ln(zT/z0)  =  1.458aH  k(S5a)'A  Pr:/'  K(vf. 

And  matching  the  slopes  there,  too, 

'A  =  1.458aH  k(85a)'/>  PrVl  re/3)'1  1(^/2)“  exp(-$s/2) 

+  i.44^[>e/,,iis)-re/,)  +  ^P(-k))\~ 


(50) 


(51) 


(52) 


(53) 


(54) 


I  again  solve  eq  54  for  ?}s  by  Newton’s  method  and  then  find  Zj/zo  from  eq  53, 

zT/zo  =  exp[-l ,458aH k(%5a)'A  Prv’  K(i h)l  (55) 

The  constant  a  is  yet  to  be  specified.  The  viscous  sublayer  velocity  profile,  eq  33,  is  approx¬ 
imately  valid  from  the  surface  to  the  lower  boundary  of  the  inertial  sublayer  at  30  v/u,  (Ten¬ 
nekes  and  Lumley  1972,  p.  160).  Therefore,  integrating  this  profile  from  zero  to  30  v/u, 
should  yield  an  average  velocity  0  for  the  viscous  sublayer.  That  average  is  0  =  15 Com¬ 
paring  this  result  with  eq  40,  we  see  that  a  should  be  roughly  15.1  have  found  that  the  value  a 
=  10  fits  the  available  data  best. 

Figure  5  shows  model  calculations  of  the  nondimensional  matching  height  as  a  function  of 


Figure  5.  Nondimensional  matching 
height  over  an  aerodynamically 
smooth  surface  as  a  function  of  a. 
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Figure  6.  Matching  of  interfacial  and  in -  Figure  7.  Current  model  predictions  for  an  aerody- 
ertial  sublayer  profiles  of  temperature  (a  namically  smooth  surface  compared  with  mea- 
=  0.71)  and  water  vapor  fa  =  0.63)  over  sured  scalar  roughness  lengths  for  water  vapor, 
an  aerodynamically  smooth  surface.  thorium-B  and  heat,  and  with  models  by  Brutsaert 

(1975b)  and  von  Karman  ( Goldstein  1965). 


o.  For  temperature  and  water  vapor,  the  value  h  u,/v  is  about  32.  Brutsaert  (1975a)  did  his 
matching  by  assuming  that  h  u,/v  =  30  for  all  values  of  a.  Figure  6  shows  the  matching  of 
interfacial  and  inertial  sublayer  profiles  for  temperature  and  water  vapor  over  an  aerody¬ 
namically  smooth  surface.  As  with  the  aerodynamically  rough  surface,  both  profiles  and 
their  first  derivatives  are  continuous  at  h. 

Finally,  Figure  7  compares  my  model  predictions  with  the  scanty  data  available  from  flows 
over  smooth  surfaces.  In  the  figure  the  data  from  Chamberlain  (1968)  are  the  “smooth  sur¬ 
face”  values  from  his  Tables  2  and  3.  The  data  from  Dipprey  and  Sabersky  (1963)  are  from 
their  smooth  pipe  (E-3,  their  Fig.  5)  and  from  their  three  rough  pipes  (their  Fig.  6-8)  for  runs 
when  R,  <  0.135.  The  fourth-order  polynomial 

ln(zs/2o)  =  0.0399  -  3.92  In*  -  1 ,22(lna)!  -  0.254(lno)J  -  0.0748(lna)4  (56) 


is  a  good  representation  of  my  model  results  for  0.35  <  a  <  10.0.  The  figure  also  shows 
Brutsaert’s  (1975b)  prediction, 


Zs/Zo  =  expHt(13.6a*  -  13.5)1,  (57) 

and  von  Karman’s  (Goldstein  1965,  p.  657;  Monin  and  Yaglom  1971,  p.  342), 

Zs/Z0  =  exp[-5  Arj(a-l)  +  ln[l  +0.83(a-l])|].  (58) 

The  three  models  are  so  close  that,  with  the  scatter  in  the  experimental  data  and  their  spars¬ 
ity,  it  is  impossible  to  decide  which  is  best. 

Notice  in  Figure  7  that  all  three  models  predict  zs/zo  =  1  for  a  =  1.  This  is  compatible 
with  the  Reynolds  analogy — that  over  an  aerodynamically  smooth  surface,  where  the  rough¬ 
ness  elements  cannot  cause  momentum  transfer  through  pressure  forces,  the  transfer  must  be 
identical  for  momentum  and  for  a  scalar  contaminant  with  a  =  1 .  von  Karman  (Goldstein 
1965,  p.  657)  explicitly  assumed  the  validity  of  the  Reynolds  analogy  and  thus  forced  his 
model  to  predict  zs/zo  =  1  at  a  =  1.  While  neither  Brutsaert  (1975a,  b)  nor  I  made  this  as¬ 
sumption,  his  model  predicts  zs/Zo  =  0.96  at  a  =  1,  and  mine  predicts  zs/zo  =  104. 


SCALAR  TRANSFER  COEFFICIENTS 

With  the  results  of  the  last  two  sections  we  can  specify  Zy/Zo  and  Zq/zo  over  snow  or  sea  ice 
for  all  /?,  (Fig.  8).  Since  temperatures  will  be  0°C  or  less,  I  take  the  Prandtl  number  (v/D)  as 
0.71  and  the  Schmidt  number  (i >/Z?w)  as  0.63,  with  values  for  the  molecular  diffusivity  of 
water  vapor,  Z>w,  taken  from  Pruppacher  and  Klett  (1978,  p.  413).  In  R .  I  evaluate  u  at  -5°C. 
Because  the  model  predicts  z%/Zo  only  over  aerodynamically  smooth  and  rough  surfaces,  to 
obtain  zs/z0  values  in  the  transition  region,  1  did  a  log-log  interpolation  between  model  re¬ 
sults  at  R,  =  0.135  and  R ,  =  2.5. 


Figure  8.  Model  predictions  of  zT/z  0  and 
Zq/zq  over  snow  and  sea  ice. 


Table  1.  Values  of  the  coefficients  in  the  polynomials 
(eq  59)  that  predict  zs/z0  for  temperature  (a  =  0.71) 
and  water  vapor  ( a  =  0.63). 

R  <  0. 135  0.135  <R  <2.5  2.5  <R  ^1000 


Temperature 


Water  vapor 


Figure  8  shows  that  Zq  is  slightly  larger  than  zT  at  all  roughness  Reynolds  numbers.  As  Fig¬ 
ures  3,  4  and  7  imply,  this  is  strictly  an  effect  of  the  difference  in  molecular  diffusivities. 
Both  zx  and  Zq  are  usually  less  than  z<j;  Zj/Zq  and  Zq/Zq  become  less  than  one  in  the  transition 
region  and  decrease  monotonically  in  the  aerodynamically  rough  region.  Thus,  zT  and  zy  are 
virtually  always  less  than  zo  in  natural  flows. 

Several  recent  models  predicted  scalar  transfer  over  aerodynamically  rough  surfaces.  The 
model  by  Garratt  and  Hicks  (1973) — which  is  just  the  empirical  equation  that  Owen  and 
Thomson  (1963)  derived — predicts  zs/Zo  values  generally  five  times  larger  than  my  model. 
The  predictions  of  zs/zo  by  Liu  et  al.  (1979),  which  admittedly  are  for  a  water  surface,  are  al¬ 
most  an  order  of  magnitude  smaller  than  mine.  Brutsaert’s  (1975b)  model  predicts  zs/zo  val¬ 
ues  smaller  than  mine  for  R ,  <  20  but  values  fairly  close  to  mine  for  larger  /?,. 

For  facilitating  computer  modeling,  I  have  fitted  the  model  results  in  Figure  8  with  poly¬ 
nomials  of  the  form 

ln(zs/zo)  =  *o  +  *i  lnfl,  +  &2(ln  /?.)*•  (59) 

Table  1  lists  the  coefficients  for  smooth,  transition  and  rough  surfaces. 

As  eq  1 1  and  12  show,  if  we  know  zT/zo,  Zq/zo  and  CD,  we  can  find  CH  and  Ct  .  The  model 
for  CD  presented  by  Banke  et  al.  (1980)  synthesizes  a  representative  collection  of  CD  values 
measured  over  a  host  of  sea  ice  surfaces.  Their  empirical  result, 

10’  CD  =  1.10  +  0.072  £,  (60) 

parameterizes  the  drag  coefficient  in  terms  of  the  rms  surface  roughness  £  in  centimetres. 
Banke  et  al.  (1980)  found  £  by  measuring  the  surface  elevation  at  1-m  intervals  for  several 
hundred  metres  upwind  of  their  instruments.  Integrating  the  power  spectrum  of  these  data 
over  wavelengths  less  than  13  m  yielded  £!.  The  zo  values  reported  by  Schmidt  (1982)  for 
blowing  snow  and  those  summarized  by  Chamberlain  (1983)  for  drifting  snow  and  sand  are 
generally  in  the  range  reported  by  Banke  et  al.  (1980);  hence,  we  assume  that  eq  60  is  a  valid 
model  for  snowfields,  too. 

Kind  (1976)  and  Chamberlain  (1983)  suggested  that  roughness  lengths  for  snow  and  sand 
obey  Charnock’s  (1955)  relation, 

ZO  =  0«J/S,  (61) 
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Figure  9.  Bulk  transfer  coefficients  for  sensible  (CH) 
and  latent  (C^j  heal  over  snow  or  sea  ice  as  a  func¬ 
tion  of  the  rms  surface  roughness  (in  centimetres ) 
and  the  IO-m  wind  speed  Uio.  The  arrows  on  the  right 
show  C D  for  the  indicated  £  value. 


where  g  is  the  acceleration  of  gravity  and  ti  is  a  dimensionless  constant  in  the  range  0.010- 
0.016.  Through  eq  10,  eq  61  would  yield  CD;  but  because  eq  60  parameterizes  the  form  drag, 
which  as  I  explained  is  an  important  effect  over  sea  ice,  I  prefer  eq  60  to  61. 

The  routine  for  estimating  CH  and  CF  is  first  to  select  a  £  value;  this  value  then  defines 
from  eq  60.  CD,  in  turn,  has  a  one-to-one  relationship  with  through  eq  10.  Finally,  we 
compute  R,  from 

R.  =  £/,.  Cg  z0/v,  (62) 

where  Ul0  is  the  wind  speed  at  10  m.  Substituting  /?,  into  eq  59,  we  use  the  resulting  zT/;o  and 
Zq/zo  values  to  compute  CH  and  CE  for  a  10-m  reference  height  from  eq  1 1  and  12.  Figure  9 
shows  CH  and  C(  as  functions  of  £  and  Ut».  Remember  again,  CD,  CH  and  C(  are  the  values 
at  neutral  stability. 


According  to  Figure  9,  CE  is  always  1-3%  larger  than  CH.  And  except  for  very  low  wind 
speed,  when  the  surface  is  aerodynamically  smooth  or  in  transition,  both  are  smaller  than 
C0. 

The  form  of  the  CH  and  CE  functions  in  Figure  9  is  different  from  that  predicted  by  Kondo 
(1975)  and  by  Liu  et  al.  (1979)  for  CH  and  CE  over  the  ocean.  According  to  Figure  9,  CH  and 
CE  are  monotonically  decreasing  functions  of  wind  speed,  except  at  very  low  speeds  over 
smooth  surfaces  where  they  are  constant.  Kondo  (1975)  predicted  that  over  the  ocean  CH  and 
CE  have  minimums  at  about  2  m/s  and  then  increase  gradually  as  the  wind  speed  increases. 
Liu  et  al.  (1979)  predicted  that  CH  and  CE  have  local  minimums  at  roughly  2  m/s,  local  maxi- 
mums  at  about  5  m/s,  then  decrease  slowly  for  increasing  wind  speeds.  The  basic  reason  for 
the  differences  between  my  model  and  these  is  that  over  the  ocean  CD  has  a  wind  speed  de¬ 
pendence.  No  such  wind  speed  dependence  has  been  established  for  the  drag  coefficient  over 
snow  or  sea  ice. 

In  Figure  9,  CH  and  CE  are  generally  between  1.0x10'*  and  1.5x10  *.  Only  over  the 
roughest  surfaces — and  then  only  at  low  wind  speeds— are  CH  and  CE  larger  than  1 .5  x  10  *. 
CH  and  CE  go  below  1.0x10'*  only  at  unusually  high  wind  speeds.  The  model  predictions 
therefore  cast  doubt  on  the  scalar  transfer  coefficients  reported  by  Hicks  and  Martin  (1972) 
and  Thorpe  et  al.  (1973).  For  wind  speeds  of  about  3  m/s,  Hicks  and  Martin  (1972)  found 
average  values  of  CH  and  CE  over  snow-covered  Lake  Mendota  to  be  0.9  x  10  *  and  2.5  x  10  ’ 
respectively.  Over  ice  in  the  Beaufort  Sea,  Thorpe  et  al.  (1973)  found  averages  of  CH  = 
1.2  x  10  ’  and  CE  =  0.55  x  10'*  for  winds  ranging  from  5  to  10  m/s.  Only  the  CH  measure¬ 
ment  by  Thorpe  et  al.  (1973)  is  compatible  with  theoretical  predictions. 

I  do  not  mean  here  to  deprecate  the  efforts  of  these  scientists;  measuring  CH  and  CE  over 
natural  snow  and  sea  ice  surfaces  is  an  important  but  extremely  difficult  job.  First,  you  have 
to  measure  «.  either  by  measuring  the  vertical  velocity  profile  or  by  measuring  t  directly. 
Next,  you  must  measure  t,  and  q, — again,  either  by  measuring  the  vertical  profiles  of  tem¬ 
perature  and  water  vapor  or  by  measuring  Hs  and  HL  directly.  These  are  necessary  not  only 
for  finding  CH  and  CE  but  also  for  making  stability  corrections.  Last,  and  probably  most  im¬ 
portant,  is  the  measurement  of  T0  and  Qq.  Since  Tq-Tt  and  Qo~Qr  are  rarely  large  over  fro¬ 
zen  surfaces,  the  T0  and  Qb  measurements  must  be  precise.  But  because  the  surface  is  ill- 
defined,  simply  deciding  what  level  over  the  snow  corresponds  to  T0  and  Q0  is  a  problem; 
finding  instruments  capable  of  measuring  T0  and  Qq  without  disturbing  the  integrity  of  the 
surface  is  another.  Consequently,  the  most  careful  flux  measurements  are  of  little  value  for 
specifying  CH  and  CE  if  T0  and  Q0  are  not  measured  as  carefully. 


CONCLUSIONS 

1  have  modeled  the  transfer  of  the  passive  scalar  contaminants  temperature  and  water 
vapor  over  aerodynamically  rough  and  smooth  snow  and  sea  ice  surfaces.  The  basis  of  the 
model  is  a  smooth  matching  of  interfacial  and  inertial  sublayer  profiles.  The  inertial  sublayer 
profile  has  the  usual  semi-logarithmic  form;  I  derive  the  interfacial  sublayer  profiles  over 
smooth  and  rough  surfaces  on  the  basis  of  a  turbulent  surface-renewal  model.  This,  I  think, 
is  the  first  such  derivation  of  the  interfacial  sublayer  profile  for  a  passive  scalar  over  an  aero¬ 
dynamically  smooth  surface. 

The  model  yields  values  of  zx/zo  and  Zy/zp  as  functions  of  the  roughness  Reynolds  num¬ 
ber.  Using  these  values  and  the  empirical  model  for  the  drag  coefficient  over  sea  ice  given  by 
Banke  et  al.  (1980),  I  predict  the  bulk  coefficients  for  sensible  (CH)  and  latent  (CE)  heat  at 
neutral  stability  over  snow  and  sea  ice.  These  depend  on  the  wind  speed  and  on  a  surface 
roughness  parameter.  CE  is  1-3%  larger  than  CH;  at  wind  speeds  greater  than  3  m/s  both  are 
virtually  always  between  1 .0  x  10" ’  and  1 .5  x  10  *.  Only  at  low  wind  speeds— which  usually  do 
not  persist— and  over  very  rough  surfaces  are  CH  and  CE  larger  than  1.5  x  10*. 
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